function [mdo, xgdo] = mpsets2(mpsd, r1, r2, t, trajdata, xgd)
% function mpsdo = mpsets2(mpsd, r1, r2, t, trajdata, xgd)
% Sets up information for each trajectory in the second stage run
% This file sets up information for each one of the individuals second stage runs
% This includes:
% 1- set up case conditions for the case (net load, committed gens, and last period schedule)
% 2- Apply additional restriction to conventional gens only induced by previous period dispatch and ramp limits
% 3- Build mpsd output struct using loadmpsd function to use as input of mpsopfl
% 
% Implementation questions and concerns:
% - xgd is being passed as input directly. Is it necessary to set InitialPg
%   to the dispatch of the previous (second stage) time period?
%   Also, de we need to set CommitKey or any of the UC parameters to any
%   particular value?
% - Presumably, xgd could be constructed directly from the input mpsd
%   since mpsd was itself built with loadmpsd using the same xgd here taken
%   as input in the tr_st function. However, this would require additional
%   coding that doesn't seem to be useful but for this specfic function.
% 
%
% Inputs
%     MPSD      :     Casefile used in first stage
%     R1        :     Structure with first stage results
%     R2        :     Structure with second stage results
%     T         :     Time period considered
%     TRAJDATA  :     Profiles struct containing modifications specific for a single trajectory to be applied to case file for period t
%     XGD       :     xGenData structure as generated by function loadxgendata
%
% Outputs
%     MPCO      :     Casefile modified with restrictions for second stage
%
% Assumptions
% (ig = conventional generators, il = demands, iw = wind units)
% - first stage results are provided
% - CommitSched are saved as a field in mpc, but do not play a role (check this)
% - The offers are taken from mpsd.offer(t)
% - Loads and stochastic generators are treated differently to conventional generators
% - Loads assumed to start no shedding (t = 1), set as mpco.gen(il, PMIN) after profile scaling
% - Full mpsd is passed (maybe just pass the mpsd for the specific time period considered)
% - Fields with information of ESS/Wind units are available mpsd.mpc.iwind, mpsd.mpc.iess (even if empty)
% - dispatches initialized in first period as r1.InitialPg
% - contingency number is set to 0 in apply_contingency - as incoming label is 0
% - contab table is taken from mpsd.cont(t, js).contab, js =1
% - Code assumes second stage time periods are equal to first stage periods
% - Running the problem, there were some nonconvergence issues. Therefore, a check to increase the range was included
% - Possibility to change energy.baseGmin to deal with convergence issues (test mpsets2_proposal1, commented)
% - Warning: This function modifies the mpc struct contained in stage 1
%   results r1 by restricting the PMAX and PMIN of conventional generators
%   (all but loads and wind gens) due to the ramp limits and the dispatch
%   at the prior time period t-1. Subsequently, a given trajectory profile
%   struct trajdata with info for a single time period is applied. Note
%   then that if the profiles where to modify the original PMIN or PMAX of
%   conventional gens, they would modify the updated values and not the
%   original PMIN and PMAX values.
%
% Verify
% When determining the ranges for Gmax and Gmin, there is an implicit constraint that may limit the transition from base(t-1) to contingency(t)
% 2013.10.23: two ranges used to deal with this: energy.Gmax/Gmin and energy.baseGmax/Gmin for base/contingency cases
%
% Future Enhancements
% - pass information on what is the previous step (s1 or internal s2) 
% - pass the contab table to be used
% - improve the changes across time (take results from s1, vs s2 in receding horizon)
% - Pass additional information for intra second period changes, e.g. 
% TRAJDATAI :     Raw data, to create information on limits for loads
% 
% 
% 2013.10.01
% Daniel Munoz
% Alberto J. Lamadrid
% 
% Modifications:
% 2014.12.12, Daniel: Second stage dispatch intervals for base cases and contingency cases are set independently.    
%                     energy.Gmax/Gmin and energy.baseGmax/Gmin are used to set those intervals separately.      
%                     Interval for base cases is set with max and min dispatches over base case scenarios.
%                     Interval for contingent cases is set using Pc, down reserves and up reserves from s1 for the specified time period. 
% 2015.05.22
% - correct  to take into account that second stage is ran using mpsopfl, not sopf2 
%   - includes location of variables such as Pactual_prior

%   MOST Paper Simulations
%   Copyright (c) 2013-2018 by Daniel Munoz-Alvarez, Alberto J. Lamadrid
%
%   This file is part of MOST Paper Simulations.
%   Covered by the 3-clause BSD License (see LICENSE file for details).

%%----- set basic constants -----
define_constants;
[CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, ...
    CT_TAREABUS, CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, ...
    CT_REP, CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, ...
    CT_LOAD_ALL_PQ, CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, ...
    CT_LOAD_FIX_P, CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, ...
    CT_MODCOST_F, CT_MODCOST_X] = idx_ct;

%%----- set basic info including commitment -----
mpco = r1.mpc;
if r1.UC.CommitKey
  mpco.gen(:, GEN_STATUS) = r1.UC.CommitSched(:, t);
end
ng = size(mpco.gen, 1);
nt = mpsd.idx.nt;
js = 1;                             % scenario from which contab table will be taken

% Obtain number of base case scenarios
nj = size(r1.flow,2);

%%----- set offers for run ----- Not needed if xgd is passed
% roffer = [
%   mpsd.offer(t).PositiveActiveReservePrice,...
%   mpsd.offer(t).PositiveActiveReserveQuantity,...
%   mpsd.offer(t).NegativeActiveReservePrice,...
%   mpsd.offer(t).NegativeActiveReserveQuantity,...
%   mpsd.offer(t).PositiveActiveDeltaPrice,...
%   mpsd.offer(t).NegativeActiveDeltaPrice
% ];
% mpco = offer2mpc(mpco, roffer);



%%----- create constraints for generators and demands -----
% Compute s2 generation ranges for all generators except wind and
% dispatchable loads (ess units not yet supported).

if t > 1
  sc = 1;
%  Pactual_prior = r2.results{t-1, sc}.gen(:, PG);
%  Pactual_prior = r2.results{t-1, sc}.results.Pc;
  Pactual_prior = r2.results{t-1, sc}.flow(1, sc, 1).mpc.gen(:, PG);
else
  Pactual_prior = r1.InitialPg;
end
RampLimit = r1.flow(t,1,1).mpc.gen(:, RAMP_30)*2*r1.Delta_T;

ig = ~isload(mpco.gen) & ~iswind(mpco.gen, mpco.iwind);

mpco.gen(ig, PMAX)   = min(mpco.gen(ig, PMAX),...
                               Pactual_prior(ig) + RampLimit(ig));
                           
mpco.gen(ig, PMIN)   = max(mpco.gen(ig, PMIN),...
                               Pactual_prior(ig) - RampLimit(ig));

if any(mpco.gen(ig, PMAX)-mpco.gen(ig, PMIN)< 0)
    unfeas = find(mpco.gen(ig, PMAX)-mpco.gen(ig, PMIN) < 0);
    if any(mpco.gen(unfeas, PMAX)-mpco.gen(unfeas, PMIN) < -0.1)
        mpco.gen(unfeas, PMIN) = mpco.gen(unfeas, PMAX);
        warning('%i significantly-unfeasible second-stage (conventional) generation range(s) was (were) modified to be feasible.',size(unfeas,1))
    else
        mpco.gen(unfeas, PMIN) = mpco.gen(unfeas, PMAX);
        warning('%i acceptably-unfeasible second-stage (conventional) generation range(s)  was (were) modified to be feasible.',size(unfeas,1));
    end
end


%%----- set initial dispatches -----
mpco.gen(:, PG) = Pactual_prior;        % Apparently, this is of no need because it does not constrain the problem at all 
il = isload(mpco.gen);
mpco.gen(il, PG) = mpco.gen(il, PMIN);

%%----- contingencies -----
contab = mpsd.cont(t, js).contab;

%%----- xGenData -----
% xgd is being passed as input
% Need to set InitialPg to Pactual_prior
xgd.InitialPg = Pactual_prior;

%%----- storage ----- Not supported yet
sd = [];

%%----- build transition prob matrix for single period single scenario -----
transmat{1} = 1;

%%----- Truncate given profiles for specified time period t
for m = 1:size(trajdata, 1)         % cycle over each profile element
  trajdata(m, 1).values = trajdata(m).values(t, :, :); % modify each profile element by selecting only current time period t
end

xgdo = xgd;
xgdo.CommitSched = mpco.gen(:, GEN_STATUS);
%xgdo.CommitKey = [];               % Creates error in loadmpsd
if isfield(xgdo, 'CommitKey')
  xgdo = rmfield(xgdo, 'CommitKey');
end
%mdo = loadmpsd(transmat, mpco, xgdo, sd, contab, trajdata);
mdo = loadmd(mpco, transmat, xgdo, sd, contab, trajdata);
mdo.UC.CommitKey = [];            % second stage does not run UC